MODULE rrf_module

#ifdef ARPREC
USE mpmodule
#endif
IMPLICIT NONE

PRIVATE

!  Machine-dependent values

!  LONG integers
!  Use "make precision COMPILER=compiler" to find the maximum integer
!  size (and real precision) available.
!  18 decimal digits is the maximum for INTEGER*8
!  gfortran's INTEGER*16 allows up to 38 digits
#ifdef GFORTRAN
INTEGER, PARAMETER, PUBLIC :: long=SELECTED_INT_KIND(38)
#else
INTEGER, PARAMETER, PUBLIC :: long=SELECTED_INT_KIND(18)
#endif
INTEGER(KIND=long), PARAMETER, PUBLIC :: max_long_int=huge(1_long)

#ifdef ARPREC
TYPE(mp_integer), PUBLIC :: one, zero
#else
INTEGER(long), PARAMETER, PUBLIC :: one=1_long, zero=0_long
#endif

INTEGER, PARAMETER :: qp=SELECTED_REAL_KIND(18)
REAL(qp) :: max_real=huge(1.0_qp)

!  Characters
CHARACTER(LEN=1), PARAMETER ::                                         &
    arrow="^", space=" ", dot=".", plus="+", minus="-", uci="I",       &
    lci="i", slash="/", star="*", shriek="!", bra="(", ket=")"

CHARACTER(LEN=1), PARAMETER :: digit(0:10)=                            &
    (/'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '0'/)

!  LIN and LOUT are the units used for READ and WRITE respectively.
INTEGER :: lin=5, lout=6

LOGICAL :: debug=.false.
LOGICAL, PUBLIC :: hard_errors=.true., error=.false.

INTEGER, PUBLIC :: nprime
INTEGER(long), ALLOCATABLE, PUBLIC :: prime(:)

TYPE node
  !  Contains the factor prime^x
  INTEGER :: x
  INTEGER(long) :: prime
  TYPE(node), POINTER :: next=>null()
END TYPE node
TYPE rrf
  TYPE(node), POINTER :: list=>null()
#ifdef ARPREC
  TYPE(mp_integer) :: factor=zero
#else
  INTEGER(long) :: factor=zero
#endif
END TYPE rrf

!  Can be set for printing in the event of an error
CHARACTER(LEN=80), PUBLIC :: args="" 

!  free points to a linked list of free nodes, or is null
TYPE(node), POINTER, SAVE :: free=>null()

PUBLIC rrf, clear, copy, rename, unit, setzero, iszero, nonzero, exch,  &
    chsign, conjugate, list, stri, init_rrf, report,                    &
    mult, power, add, subtr, root, int_to_rrf, multi, normalize, equal, &
    pp_to_rrf, number, display, rrf_to_real, rrf_to_4i, char4i,         &
    write4i, rrf_to_char, exp_to_char, setf, multf

INTERFACE int_to_rrf
#ifdef ARPREC
  module procedure fromi_int, fromi_mp
#else
  module procedure fromi_int, fromi_long
#endif
END INTERFACE int_to_rrf

INTERFACE number
  MODULE PROCEDURE number4, number_long
END INTERFACE

!  Factorial table
INTEGER, PARAMETER, PUBLIC :: MAXFACT=1000
TYPE fctorial
  LOGICAL :: unset=.true.
  TYPE(rrf), POINTER :: rrf=>null()
END type fctorial
TYPE(fctorial), SAVE, PUBLIC :: factorial(0:MAXFACT)

CHARACTER(LEN=*), PARAMETER :: this_file="rrf_module.F90"

CONTAINS

SUBROUTINE init_rrf(quiet,primes)

LOGICAL, INTENT(IN), OPTIONAL :: quiet
INTEGER, INTENT(IN), OPTIONAL :: primes

INTEGER(long) :: m
INTEGER :: i, j, ok
LOGICAL :: verbose=.true.

#ifdef ARPREC
one=1
zero=0
#endif

!  Initialise factorial table (initially contains just 0! and 1!).
call unit(factorial(0)%rrf)
call unit(factorial(1)%rrf)
factorial(0)%unset=.false.
factorial(1)%unset=.false.


if (present(quiet)) then
  verbose=.not. quiet
else
  verbose=.true.
end if
if (present(primes)) then 
  nprime=primes
else
  nprime=200000
endif

if (verbose) then
  print "(a/a//a)", "Root rational fraction program",                  &
      "Version 4.0.3  May 2013"
  print "(a)", "Constructing table of primes ..."
end if

!  Set up table of first NPRIME primes
allocate(prime(nprime),stat=ok)
if (ok>0) then
  print "(a,i0,a)", "Can't allocate ", nprime, " primes"
end if
prime(1)=2
m=1
j=1
do
  m=m+2
  do i=1,j
    if (mod(m,prime(i)) == 0) exit   ! not prime
    if (m < prime(i)*prime(i)) then  ! prime found
      j=j+1
      prime(j)=m
      exit
    endif
  end do
  if (j == nprime) exit
end do

END SUBROUTINE init_rrf

!-----------------------------------------------------------------------

SUBROUTINE get(p)
!  Get a single link and point p at it

TYPE(node), POINTER :: p
INTEGER :: ok

if (associated(free)) then
  p=>free
  free=>free%next
else
  allocate(p,stat=ok)
  if (ok>0) print "(a)", "Couldn't allocate node"
end if

END SUBROUTINE get

!-----------------------------------------------------------------------

SUBROUTINE link(p)
!  Insert a new link after p and point p at it

TYPE(node), POINTER :: p
TYPE(node), POINTER :: pz=>null()

call get(pz)
pz%next=>p%next
p%next=>pz
p=>pz
pz=>null()

END SUBROUTINE link

!-----------------------------------------------------------------------

SUBROUTINE unlink(p)
!  Remove link p%next and point p at the one after that

TYPE(node), POINTER :: p
TYPE(node), POINTER :: pz=>null()

pz=>p%next
p%next=>pz%next
pz%next=>free
free=>pz
pz=>null()

END SUBROUTINE unlink

!-----------------------------------------------------------------------

SUBROUTINE clear(k)
!  If rrf k doesn't exist, allocate it.
!  Otherwise free k%list and set k%factor=0

TYPE(rrf), POINTER :: k
TYPE(node), POINTER :: i=>null(), j=>null()

if (.not. associated(k)) then
  allocate(k)
else
  if (associated(k%list)) then
    j=>k%list
    if (associated(free)) then
      i=>free
      free=>j%next
      j%next=>i
      i=>null()
    else
      free=>j%next
      j%next=>null()
    end if
  endif
  k%list=>null()
  k%factor=zero
  j=>null()
end if

END SUBROUTINE clear

!-----------------------------------------------------------------------

SUBROUTINE copy(k,kz)
!  Copy rrf K to KZ

TYPE(rrf), POINTER :: k, kz
TYPE(node), POINTER :: p, pz

if (.not. associated(k)) call report("Source for COPY is undefined")
if (associated(k,kz)) call report("Source and target of COPY are identical")

call clear(kz)
if (iszero(k)) then
  call setzero(kz)
else
  call unit(kz)
  pz=>kz%list
  p=>k%list
  pz%x=p%x
  do
    p=>p%next
    if (p%prime < 0) exit
    call link(pz)
    pz%prime=p%prime
    pz%x=p%x
  end do
  kz%factor=k%factor
end if

END SUBROUTINE copy

!-----------------------------------------------------------------------

SUBROUTINE rename(k1,k2)

TYPE(rrf), POINTER :: k1, k2

call clear(k2)
k2=>k1
k1=>null()

END SUBROUTINE rename

!-----------------------------------------------------------------------

SUBROUTINE unit(k)
!  Set K to unity. Described by an rrf with factor 1 and the single
!  list entry [(-1)^0] -> itself


TYPE(rrf), POINTER :: k
TYPE(node), POINTER :: p=>null()

! print "(a)", "In unit subroutine"
call clear(k)
! print "(a)", "Cleared"
call get(p)
k%list=>p
p%next=>p
p%prime=-1
p%x=0
k%factor=one

END SUBROUTINE unit

!-----------------------------------------------------------------------

SUBROUTINE setzero(k)
!  Set K to zero. Described by an rrf with factor 0 and the single
!  list entry [(-1)^0] -> itself

TYPE(rrf), POINTER :: k
TYPE(node), POINTER :: p

call clear(k)
call get(p)
k%list=>p
p%next=>p
p%prime=-1
p%x=0
k%factor=zero

END SUBROUTINE setzero

!-----------------------------------------------------------------------

LOGICAL FUNCTION iszero(k)

TYPE(rrf), POINTER :: k

if (.not. associated(k)) then
  call report("Undefined argument to iszero")
else
  iszero=(k%factor==zero)
end if

END FUNCTION iszero

!-----------------------------------------------------------------------

LOGICAL FUNCTION nonzero(k)

TYPE(rrf), POINTER :: k

if (.not. associated(k)) then
  call report("Undefined argument to nonzero")
else
  nonzero=(k%factor.ne.zero)
end if

END FUNCTION nonzero

!-----------------------------------------------------------------------

SUBROUTINE exch(k1,k2)

TYPE(rrf), POINTER :: k1, k2
TYPE(rrf), POINTER :: i=>null()
i=>k1
k1=>k2
k2=>i
i=>null()
END SUBROUTINE exch

!-----------------------------------------------------------------------

SUBROUTINE chsign(k)

TYPE(rrf), POINTER :: k
TYPE(node), POINTER :: p=>null()

if (k%factor == zero) return
p=>k%list
p%x=modulo(p%x+2,4)
p=>null()

END SUBROUTINE chsign

!-----------------------------------------------------------------------

SUBROUTINE conjugate(k)

TYPE(rrf), POINTER :: k
TYPE(node), POINTER :: p=>null()

if (k%factor == zero) return
p=>k%list
p%x=modulo(4-p%x,4)
p=>null()

END SUBROUTINE conjugate

!-----------------------------------------------------------------------

SUBROUTINE list(k,string)
!  Print rrf K

TYPE(rrf), POINTER :: k
CHARACTER(LEN=*), OPTIONAL :: string

TYPE(node), POINTER :: p=>null()

if (present(string)) then
  print "(a)", string
end if
if (.not. associated(k)) then
  print "(a)", '  undefined'
else if (k%factor==zero) then
  print "(a)", '  zero'
else
  p=>k%list
  do
    print "(i18,i5)", p%prime, p%x
    p=>p%next
    if (associated(p)) then
      if (associated(p,k%list)) then
        !  Circular list -- back to start
        !  print "(1x)"
        exit
      endif
    else
      print "(a)", "         => null"
      exit
    end if
  end do
!  if (k%factor .ne. one) then
    write (6,"(16x,a)",advance="no") "* "
#ifdef ARPREC
    call mpwrite(6,k%factor)
#else
    print "(i0)", k%factor
#endif
  end if
!  end if

END SUBROUTINE list

!-----------------------------------------------------------------------

SUBROUTINE report(message,line,file)

!  Print message m; then stop if hard_errors or error is true.
!  Otherwise set error true.

CHARACTER(LEN=*), INTENT(IN) :: message
INTEGER, OPTIONAL :: line
CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: file

if (args .ne. "") write(lout,"(a)") trim(args)
if (present(line)) then
  write (lout,"(/2a,i0,2a)") trim(message), " at line ", line, " in ", trim(file)
else
  write (lout,"(/a)") message
end if

if (hard_errors .or. error) then
  write (lout,'(/1x,a)')                                                  &
      '**********       Program terminated       **********'
  stop
else
  error=.true.
endif

END SUBROUTINE report

!-----------------------------------------------------------------------

FUNCTION stri(i)
!  Return a left-justified character string representation of the
!  integer I

INTEGER, INTENT(IN) :: i
CHARACTER(LEN=20) :: stri, buffer

write (unit=buffer,fmt="(i0)") i
stri=adjustl(buffer)

END FUNCTION stri

!-----------------------------------------------------------------------

SUBROUTINE mult(k1,k2,mexp)
!  Multiply k1 by k2**mexp. mexp is optional, default value 1.

TYPE(rrf), POINTER :: k1, k2
INTEGER, INTENT(IN), OPTIONAL :: mexp

TYPE(node), POINTER :: p1=>null(), p2=>null(), p1p=>null()
INTEGER :: mx
CHARACTER(LEN=80) :: string

if (present(mexp)) then
  mx=mexp
else
  mx=1
end if

!  Special cases
if (iszero(k1) .or. mx == 0) return
if (associated(k1,k2)) then
  string="Identical arguments to MULT"
  call report(string)
endif
if (iszero(k2)) then
  if (mx < 0) then
    call report('Division by zero')
    call clear(k1)
  else
    call setzero(k1)
  endif
  return
endif

!  General case
!  p1p is the link previous to p
p1p=>k1%list
p1=>p1p%next
p2=>k2%list%next

do
  ! print "(a,i0,a,i0)", "p1%prime = ", p1%prime, ", p2%prime = ", p2%prime
  if (p1%prime == p2%prime) then
!  Term present in both K1 and K2: form product term
    p1%x=p1%x+mx*p2%x
    if (p1%prime < 0) then
      !  This is the sign factor
      p1%x=modulo(p1%x,4)
      exit
    end if
    if (p1%x == 0) then
!  Product term has zero exponent -- remove it and return link to free list
!  This removes the link _after_ p1p, i.e. p1
      call unlink(p1p)
      p1=>p1p%next
    else
      p1p=>p1
      p1=>p1%next
    endif
    p2=>p2%next
  else if (p2%prime > 0 .and.                                           &
      (p1%prime < 0 .or. p1%prime > p2%prime)) then
!  Term in K2 not in K1: insert it before p1
    call link(p1p)
    p1=>p1p
    p1%prime=p2%prime
    p1%x=mx*p2%x
    p2=>p2%next
    p1=>p1%next
    cycle
  else if (p1%prime > 0 .and.                                           &
      (p2%prime < 0 .or. p2%prime > p1%prime)) then
!  Term in K1 not in K2: advance to next term
    p1p=>p1
    p1=>p1%next
  endif
end do

!  Deal with any unfactorized factors
k1%factor=k1%factor*k2%factor

END SUBROUTINE mult

!-----------------------------------------------------------------------

SUBROUTINE power(k, m)
!  Replace K by K**M

TYPE(rrf), POINTER :: k
INTEGER, INTENT(IN) :: m

TYPE(node), POINTER :: p=>null(), prev=>null()
INTEGER :: i
INTEGER(long) :: f

if (m == 0) then
  call unit(k)
else if (m == 1) then
  return
else if (iszero(k)) then
  if (m < 0) call report('Division by zero')
else
  p=>k%list
  do
    prev=>p
    p=>p%next
    p%x=m*p%x
    if (p%prime < 0) exit
  end do
endif

if (k%factor .ne. one) then
  if (m > 0) then
    f=k%factor
    do i=1,m
      !  Check for overflow
      if (huge(1_long)/k%factor < f) then
        k%factor=k%factor*f
      else
        call report("Integer overflow",__LINE__,this_file)
      end if
    end do
  else if (m == 0) then
    k%factor=one
  else
    !  Have to insert negative power of factor as power of "prime"
    !  This will probably lead to trouble later
    call get(prev)
    prev%prime=k%factor
    prev%x=2*m
    print "(a,i0,a)", "Warning: unfactorizable number ", k%factor,      &
        " raised to negative power"
    k%factor=one
  end if
  prev=>null()
end if

END SUBROUTINE power

!-----------------------------------------------------------------------

SUBROUTINE add(k1, k2, debug)
!  Add K2 to K1
!  Normally k1%factor will contain 1 on entry and exit.  If however it proves
!  impossible to decompose the sum into primes small enough to occupy
!  one integer location, the oversize factor is left in k1%factor.
!  At the end of a sequence of additions the result can usually be
!  factorised satisfactorily.

TYPE(rrf), POINTER :: k1, k2
LOGICAL, OPTIONAL, INTENT(IN) :: debug

LOGICAL :: next1, next2, dbg
TYPE(rrf), POINTER :: k=>null(), k10=>null(), k20=>null()
TYPE(node), POINTER :: p=>null(), p1=>null(), p2=>null()
INTEGER :: ld, m
#ifdef ARPREC
TYPE(mp_integer) :: a1, a2, pp, ppld
#else
INTEGER(long) :: a1, a2, pp, ppld
#endif
CHARACTER(LEN=200) :: buffer=""

! print "(a)", "Entering ADD routine"
if (present(debug)) then
  dbg=debug
else
  dbg=.false.
end if

!  Special cases
if (associated(k1,k2)) then
  call report('Identical arguments to ADD')
else if (iszero(k1)) then
  call copy (k2, k1)
  return
else if (iszero(k2)) then
!  K2=0:  force factorisation of K1%factor if possible
  if (k1%factor == one) return
  call normalize(k1)
  return
end if

call unit(k)
a1=k1%factor
a2=k2%factor
next1=.true.
next2=.true.
p1=>k1%list
p2=>k2%list
p=>k%list

!  The next section constructs the highest common factor of K1 and
!  K2 and stores it as K.  The multipliers are A1 and A2, so that
!  K1=A1*K and K2=A2*K.

! print *, "Starting ADD loop"
do
  if (next1) p1=>p1%next
  if (next2) p2=>p2%next
  ! print "(6(a,i0))", "p1: ", p1%prime, "^", p1%x, " a1 = ", a1,          &
  !     "  p2: ", p2%prime, "^", p2%x, " a2 = ", a2 
  if (p1%prime == p2%prime) then
    !  Entry in both K1 and K2
    if (p1%prime > 0) then
      call link(p)
      p%prime=p1%prime
    else
      p=>p%next
    end if
    p%x=min(p1%x,p2%x)
#ifdef ARPREC
    pp=mpint(p2%prime)
#else
    pp=p2%prime
#endif
    ld=p2%x-p1%x
    ! print "(4(a,i0))", "p: ", p%prime, "^", p%x,  "   ld = ", ld
    next1=.true.
    next2=.true.
  else if ((0 < p2%prime .and. p2%prime < p1%prime) .or. p1%prime < 0) then
    !  Entry in K2 only
    ! print "(4(a,i0))", "p2: ", p2%prime, "^", p2%x 
    if (p2%x<0) then
      call link(p)
      p%prime=p2%prime
      p%x=p2%x
    end if
#ifdef ARPREC
    pp=mpint(p2%prime)
#else
    pp=p2%prime
#endif
    ld=p2%x
    ! print "(4(a,i0))", "p: ", p%prime, "^", p%x,  "   ld = ", ld
    next2=.true.
    next1=.false.
  else if ((0 < p1%prime .and. p1%prime < p2%prime) .or. p2%prime < 0) then
    !  Entry in K1 only
    ! print "(4(a,i0))", "p1: ", p1%prime, "^", p1%x 
    if (p1%x<0) then
      !  Transfer this entry to the hcf
      call link(p)
      p%prime=p1%prime
      p%x=p1%x
    end if
#ifdef ARPREC
    pp=mpint(p1%prime)
#else
    pp=p1%prime
#endif
    !  Note (twice) power of p1%prime to be divided into a1
    ld=-p1%x
    ! print "(4(a,i0))", "p: ", p%prime, "^", p%x,  "   ld = ", ld
    next2=.false.
    next1=.true.
  end if
  ! print "(3i3,2f8.0)", k%prime, k%x, ld, a1, a2

  ! print "(3i6)", p%prime, p%x, ld
  if (ld .ne. 0) then
    if (mod(ld,2) .ne. 0) then
      call report('Incompatible summands')
      call list(k10)
      call list(k20)
      stop
    end if
    ld=ld/2
    !  Note: if pp<0 then pp=-1.
    ppld=pp**(abs(ld))
    if (ld > 0) then
      if (max_long_int/a2>ppld) then
        a2=a2*ppld
      else
        call overflow(a2,pp,ld)
      end if
    else
      if (max_long_int/a1>ppld) then
        a1=a1*ppld
      else
        call overflow(a1,pp,abs(ld))
      end if
    end if
  end if
  ! print "(a,i0,a,i0)", "a1 = ", a1, "  a2 = ", a2
  if (p1%prime < 0 .and. p2%prime < 0) exit
  ! call list(k)

!  End of decomposition into h.c.f. and multipliers
end do
! print *, "Decomposition done"

call setzero(k1)
!  Test for overflow
if (sign(1_long,a1) == sign(1_long,a2)) then
  if ((a2>0 .and. a1>max_long_int-a2) .or.                              &
      (a2<0 .and. a1<-max_long_int-a2)) then
    print "(a,i0,a,i0)", "Integer overflow in ADD: ", a1, " + ", a2  
    call report("Stopping")
  end if
end if
k1%factor=a1+a2
if (dbg) then
  buffer="h.c.f.:"
  m=9
  call rrf_to_char(k,buffer,m,0)
  print "(a)", trim(buffer)
  print "(3(a,i0))", "a1 = ", a1, "  a2 = ", a2, "  factor = ", k1%factor
end if
! call list(k1)
! print *, "Factorizing"
call normalize(k1,dbg)
! if (dbg) call display(k1,1)
! print *, "Factorized"
! print "(3(a,i0))", "reduced factor = ", k1%factor
! call list(k1,"k1")
! call list(k,"k")
call mult(k1,k,1)
! if (dbg) call display(k1,1)
! call list(k1,"k1")
call clear(k)

END SUBROUTINE add

!-----------------------------------------------------------------------

SUBROUTINE subtr(k1,k2)
!  Subtract K2 from K1

TYPE(rrf), POINTER :: k1, k2

call chsign(k2)
call add(k1,k2)
call chsign(k2)

END SUBROUTINE subtr

!-----------------------------------------------------------------------

SUBROUTINE overflow(a,p,x)

INTEGER(long) :: a, p
INTEGER :: x

print "(a,i0,a,i0,a,i0)", "Integer overflow in ADD: ", a, " * ", p, "^", x
call report("Stopping")

END SUBROUTINE overflow

!-----------------------------------------------------------------------

SUBROUTINE normalize(k,debug)

!  Attempt to factorize k%factor if not 1 or 0

TYPE(rrf), POINTER :: k, ka=>null()
LOGICAL, OPTIONAL :: debug
TYPE(node), POINTER :: p=>null()

!  Standardize sign representation
p=>k%list
p%x=modulo(p%x,4)
!  if abs(k%factor)>1, attempt to factorize it
if (k%factor == zero .or. k%factor == 1) return
call int_to_rrf(k%factor, ka, debug)
k%factor=one
call mult(k,ka)
if (k%factor > one) then
  if (present(debug)) then
    if (debug) then
      print '(a)', 'Warning: intermediate value not factorized: '
#ifdef ARPREC
      call mpwrite(6,k%factor)
#else
      print "(i0)", k%factor
#endif
    end if
  end if
end if

END SUBROUTINE normalize

!-----------------------------------------------------------------------

LOGICAL FUNCTION equal(k1,k2)

TYPE(rrf), POINTER :: k1, k2
TYPE(node), POINTER :: p1=>null(), p2=>null()

if (k1%factor==zero .and. k2%factor==zero) then
  equal=.true.
  !  print "(a)", "both zero"
else if (k1%factor==zero .or. k2%factor==zero) then
  !  One is null, one not
  equal=.false.
  !  print "(a,i0,1x,i0)", "one zero, one not: ", k1%factor, k2%factor
else
  call normalize(k1)
  call normalize(k2)
  p1=>k1%list
  p2=>k2%list
  do
    p1=>p1%next
    p2=>p2%next
    !  print "(i6, i3, i6, i3)", p1%prime, p1%x, p2%prime, p2%x
    if (p1%prime .ne. p2%prime .or. p1%x .ne. p2%x) then
      equal=.false.
      exit
    else if (p1%prime < 0) then
      equal=.true.
      exit
    end if
  end do
end if

p1=>null()
p2=>null()

END FUNCTION equal

!-----------------------------------------------------------------------

SUBROUTINE root(k)

TYPE(rrf), POINTER :: k
TYPE(node), POINTER :: p=>null(), prev=>null()

LOGICAL :: ok

if (iszero(k)) return
p=>k%list

p%x=modulo(p%x,4)
ok=.true.
do
  prev=>p
  p=>p%next
  if (mod(p%x,2) .eq. 0) then
    p%x=p%x/2
  else
    ok=.false.
  endif
  if (p%prime < 0) exit
end do

if (k%factor > one) then
  !  Treat factor as "prime"
  call get(prev)
  prev%prime=k%factor
  prev%x=1
  k%factor=one
  prev=>null()
end if

if (.not. ok) then
  call report('Square root argument invalid')
  ! call clear(k)
end if

END SUBROUTINE root

!-----------------------------------------------------------------------

#ifdef ARPREC
SUBROUTINE fromi_mp(n,k)
!  Express mp integer N as PP K

TYPE(mp_integer), INTENT(IN) :: n
TYPE(rrf), POINTER :: k
TYPE(node), POINTER :: q

TYPE(mp_integer) :: m
INTEGER :: p
LOGICAL :: new

call setzero(k)
if (n .eq. zero) return

if (debug) then
  print "(a)", "Factorising "
  call mpwrite(lout,n)
end if
m=n
call unit(k)
q=>k%list
if (n < zero) then
  m=-m
  q%x=2
endif
if (m == one) return

if (debug) print "(2i5,3x,i0)", q%prime, q%x, m
do p=1,nprime
  new=.true.
  do while ((m/prime(p))*prime(p) == m)
    if (new) then
      call link(q)
      q%x=0
      q%prime=prime(p)
      new=.false.
    end if
    q%x=q%x+2
    m=m/prime(p)
  end do
  if (debug .and. .not. new .and. q%x>0) print "(2i5,3x,i0)", q%prime, q%x, m
  if (m == one) exit
  if (m - prime(p)*prime(p) < zero) then
    !  Remaining factor must be prime
    call link(q)
    q%x=2
    q%prime=m
    if (debug) print "(2i5,3x,i0)", q%prime, q%x, m
    m=1
    exit
  end if
end do

if (m>one) then
  !  Treat remaining factor as unfactorized
  k%factor=m
  if (debug) then
    print "(a)", "Unfactorized: "
    call mpwrite(lout, k%factor)
end if

q=>q%next
if (debug) call list(k)

END SUBROUTINE fromi_mp

!-----------------------------------------------------------------------

#else
SUBROUTINE fromi_long(n,k,dbg)
!  Express long integer N as PP K

INTEGER(long), INTENT(IN) :: n
TYPE(rrf), POINTER :: k
LOGICAL, OPTIONAL :: dbg

TYPE(node), POINTER :: q=>null()
INTEGER(long) :: m
INTEGER :: p
LOGICAL :: new, debug=.false.

if (present(dbg) .and. abs(n) > 10000_long) then
  debug=dbg
else
  debug=.false.
end if

call setzero(k)
if (n == zero) return

if (debug) print "(a,i0)", "Factorising ", n
m=n
call unit(k)
q=>k%list
if (n < zero) then
  m=-m
  q%x=2
endif
if (m == one) return

if (debug) print "(i18,i3,3x,i0)", q%prime, q%x, m
do p=1,nprime
  new=.true.
  do while (mod(m,prime(p)) == 0)
    if (new) then
      call link(q)
      q%x=0
      q%prime=prime(p)
      new=.false.
    end if
    q%x=q%x+2
    m=m/prime(p)
  end do
  if (debug .and. .not. new .and. q%x>0) print "(i18,i3,3x,i0)", q%prime, q%x, m
  if (m == one) exit
  if (m < prime(p)*prime(p)) then
    !  Remaining factor must be prime
    call link(q)
    q%x=2
    q%prime=m
    m=1
    if (debug) print "(i18,i3,3x,i0)", q%prime, q%x, m
    exit
  end if
end do

if (m>one) then
  !  Treat remaining factor as unfactorized
  k%factor=m
  if (debug) print "(a,i0)", "Unfactorized: ", k%factor
end if

q=>q%next
if (debug) call list(k)

END SUBROUTINE fromi_long

!-----------------------------------------------------------------------

#endif
SUBROUTINE fromi_int(n, k)
!  Express integer n as RRF K.

INTEGER, INTENT(IN) :: n
TYPE(rrf), POINTER :: k

#ifdef ARPREC
TYPE(mp_integer) :: m
#else
INTEGER(long) :: m
#endif

m=n
#ifdef ARPREC
call fromi_mp(m,k)
#else
call fromi_long(m,k)
#endif

END SUBROUTINE fromi_int

!-----------------------------------------------------------------------

SUBROUTINE multi(k,i,m)
!  Multiply K by I**M, I, M integers

TYPE(rrf), POINTER :: k
INTEGER(long), INTENT(IN) :: i
INTEGER, INTENT(IN) :: m
TYPE(rrf), POINTER :: kz

kz=>null()
call fromi_long(i,kz)
call mult(k,kz,m)
call clear(kz)

END SUBROUTINE multi

!-----------------------------------------------------------------------

!  Conversion from character representations

SUBROUTINE pp_to_rrf(m,l, np, k)

!  Read an RRF from a character buffer M, starting at position L+1.
!  The number must be in power-of-prime form

!  sign  x1 x2 ...  xNP [.] n[!][^exp] [.] n[!][^exp] ...

!  where sign is +, [+]i, -, or -i, exp is [+]n[/2] or -n[/2],
!  and n is an unsigned integer.
!  In all the above expressions, brackets denote optional items and
!  do not appear in the actual input. Spaces are required as shown, but
!  two or more spaces may be replaced by a single space.
!  x1 ... xNP are the exponents of the first NP primes. (NP may be zero)
!  Subsequent integer and factorial terms may occur in any order and may
!  be interspersed with each other.
!  The ^exp may be omitted if exp is unity; the ^ may be omitted if
!  the next character is + or -.

CHARACTER(LEN=*), INTENT(IN) :: m
INTEGER, INTENT(INOUT) :: l
INTEGER, INTENT(IN) :: np
TYPE(rrf), POINTER :: k

TYPE(node), POINTER :: q=>null()
INTEGER(long) :: nn
INTEGER :: n, x, p
LOGICAL :: sign
CHARACTER :: c
TYPE(rrf), POINTER :: kz

kz=>null()
call unit(k)
q=>k%list

sign=.false.
do
  if (m(l:) == "") then
    if (sign) then
      return
    else
      call report('Empty power-of-primes expression')
      call clear(k)
      return
    endif
  end if
  l=l+1

  c=m(l:l)
  ! print "(3a)", '"', c, '"'
  select case(c)
  case(space)
!  Space
    if (sign) exit
  case(plus,minus)
    if (sign) then
      call report('Syntax error in input')
      call clear(k)
      return
    endif
    if (c == minus) then
      q%x=q%x+2
    endif
    sign=.true.
  case(uci,lci)
!  I or i
    q%x=q%x+1
    sign=.true.
!   exit
  case(digit(0))
!  0
    call clear(k)
    return
!  case(digit(1))
!  1
!    if (m(l+1:l+1) .eq. space .or. m(l+1:l+1) .eq. dot) then
!      l=l+1
!    endif
!    exit
  case default
    call report('Sign factor missing in power-of-primes expression')
    exit
  end select
end do

!  Exponents of first NP primes
do p=1,np
  x=2
  call number(m,l, nn,x)
  if (x .eq. 0) exit
  if (n .eq. 0) cycle
  call link(q)
  q%prime=prime(p)
  q%x=(nn+nn)/x
end do
q=>q%next

!  Explicit n^exponent or n!^exponent
! l=l-1
do
  if (m(l:) == "") exit
  do while (m(l:l) .eq. space .or. m(l:l) .eq. dot)
    l=l+1
  end do
  x=1
  call number (m,l, nn,x)
  if (x .eq. 0) return
  if (m(l:l) .eq. shriek) then
    n=nn
    call setf(n,kz)
    l=l+1
  else
    call int_to_rrf(nn,kz)
  endif
  select case(m(l:l))
  case(arrow,plus,minus)
    if (m(l:l) == arrow) l=l+1
    x=2
    call number(m,l, nn,x)
    if (x .eq. 0) call report('Syntax error in power-of-primes expression')
    if (x .eq. 2) call root(kz)
    n=nn
    call mult(k,kz,n)
  case(space)
    call mult(k,kz,1)
  case default
    call report('Syntax error in power-of-primes expression')
  end select
end do

call clear(kz)

END SUBROUTINE pp_to_rrf

!-----------------------------------------------------------------------

SUBROUTINE number_long(m,l, n,x)
!  If X=1 on entry, read an integer [-]nnn... terminated by any
!  non-digit;  if X has any other value on entry, read an integer or
!  half-integer [-]nnn...[/2].
!  In either case, the number is read from the buffer M, starting at
!  position L and terminated by space, non-numerical character or end of
!  buffer. L is updated to point at the next character after the number.
!  If no digit is found before a terminator, or if the digit is missing
!  after /, X is set to zero; otherwise X is set to x if /x is
!  present and expected and to 1 otherwise.

CHARACTER(LEN=*), INTENT(IN) :: m
INTEGER(long), INTENT(OUT) :: n
INTEGER, INTENT(INOUT) :: l, x

INTEGER :: i, sign
LOGICAL :: found, whole

whole=x .eq. 1
n=0
x=1
sign=1
found=.false.


if (debug) print "(a,a,i3/a,i2)", trim(m), "  length ", len(m), "Number: l=", l
l=l-1
chars: do
  l=l+1
  if (debug) print "(i2,1x,a,1x,i0)", l, m(l:l), n
  ! if (l .gt. len(trim(m)) .and. found) exit

  do i=0,9
    if (m(l:l) .eq. digit(i)) then
      !  Digit found
      if (x .eq. 0) then
        x=i
        l=l+1
        exit
      else if (n .gt. max_long_int/10) then
        call report('Integer overflow')
        x=0
        exit
      else
        found=.true.
        n=n*10+i
        cycle chars
      endif
    end if
  end do

  !  Non-digit
  select case (m(l:l))
  case(space)
    if (x .eq. 0) call report('Denominator missing')
    if (found) exit
  case(plus)
    if (found) exit
    sign=1
  case(minus)
    if (found) exit
    sign=-1
  case(slash)
    if (.not. found) then
      call report('Numerator missing')
      x=0
      exit
    else if (whole)then
      exit
    else
      x=0
      if (l .ge. len(trim(m))) then
        call report('Denominator missing')
        exit
      end if
      cycle
    end if
  case default
    if (found) then
      exit
    else
      call report('Non-digit while attempting to read number')
      x=0
      exit
    endif
  end select
end do chars

n=n*sign
if (debug) print "(a,i0)", "n = ", n

END SUBROUTINE number_long

!-----------------------------------------------------------------------

SUBROUTINE number4(m,l, n,x)

IMPLICIT NONE
CHARACTER(LEN=*), INTENT(IN) :: m
INTEGER, INTENT(OUT) :: n
INTEGER, INTENT(INOUT) :: l, x
INTEGER(long) :: nn

call number_long(m,l, nn,x)
n=nn

END SUBROUTINE number4

!-----------------------------------------------------------------------

!  Conversion from rrf to other representations

SUBROUTINE display(k,mode,label)

!  Display: mode = 1: power-of prime
!                  2: real
!                  3: (i1/i2)*sqrt(i3/i4)

TYPE(rrf), POINTER :: k
INTEGER, INTENT(IN) :: mode
CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: label

CHARACTER :: c
CHARACTER(LEN=1024) :: iout
REAL(qp) :: a
INTEGER :: j, m

! print "(a)", "Entering display"

if (.not. associated(k)) then
  call report("Undefined rrf in display")
end if

if (present(label)) then
  iout=label
  m=len_trim(label)+2
else
  m=1
  iout=""
end if

select case(mode)
case(1)
  call rrf_to_char(k, iout, m, 0)
  if (.not. error) write (lout,"(a)") iout(1:m)
  
case(2)
  call rrf_to_real(k,a,j)
  if (nonzero(k) .and. a==0.0_qp) then
    write(lout,"(8x,a)") "Real value too big"
  else
    if (j .eq. 1) then
      c = lci
    else
      c = space
    endif
    write(lout,"(g30.16,1x,a1)") a, c
  end if
  
case(3)
  call write4i(k, iout, m,132)
end select

error=.false.

END SUBROUTINE display

!-----------------------------------------------------------------------

SUBROUTINE rrf_to_real (k,a,i)
!  Express K as A*i**I, where I=0 or 1 and A is real
!  If the real is going to be too big, A is set to zero and an error
!  report is generated.

TYPE(rrf), POINTER :: k
INTEGER, INTENT(OUT) :: i
REAL(KIND=qp), INTENT(OUT) :: a

TYPE(node), POINTER :: p=>null()
REAL(KIND=qp) :: an

i=0
if (iszero(k)) then
  a=0d0
  return
endif

a=1d0

p=>k%list
do
  p=>p%next
  if (p%prime .gt. 0) then
    an=p%prime
    if (p%x > 0 .and. max_real/(an**p%x) < a) then
      call report("Real overflow")
      a=0
      return
    else
      a=a*an**p%x
    end if
  else
    a=sqrt(a)
    p%x=mod(p%x,4)
    if (p%x < 0) p%x=p%x+4
    if (mod(p%x,2) .ne. 0) i=1
    if (p%x .ge. 2) a=-a
    exit
  end if
end do

END SUBROUTINE rrf_to_real

!-----------------------------------------------------------------------

SUBROUTINE rrf_to_4i(k0,i1,i2,i3,i4)
!  Express PP K0 as (I1/I2)sqrt(I3/I4).  I2, I4 >0; I1, I3 may be
!  negative.
!  If integer overflow occurs, I2 is set to zero and I1 and I3 contain
!  rubbish.
!  Note that an attempt is made to factorize k0%factor if necessary, 
!  but if this is not completely successful the remaining factor will
!  be left in k0%factor and not accounted for in I1--I4.

TYPE(rrf), POINTER :: k0
INTEGER(long), INTENT(OUT) :: i1,i2,i3,i4

TYPE(node), POINTER :: p=>null()
INTEGER :: n
LOGICAL :: ovflow=.false.

if (k0%factor>1) then
  call normalize(k0)
end if

p=>k0%list
i1=1
i2=1
i3=1
i4=1

if (iszero(k0)) then
  i1=0
  return
end if

!  Sign
p%x=mod(p%x,4)
if (p%x < 0) p%x=p%x+4
if (mod(p%x,2) .ne. 0) i3=-1
if (p%x .ge. 2) i1=-1

ovflow=.false.
do
  p=>p%next
  if (p%prime < 0) exit
  n=p%x
  if (n .gt. 0) then
    if (mod(n,2) .ne. 0) then
      call imult(i3,int(p%prime,long),1)
      if (i3 .eq. 0) ovflow=.true.
      n=n-1
    endif
    n=n/2
    call imult (i1,int(p%prime,long),n)
    if (i1 .eq. 0) ovflow=.true.
  else
    n=-n
    if (mod(n,2) .ne. 0) then
      call imult(i4,int(p%prime,long),1)
      if (i4 .eq. 0) ovflow=.true.
      n=n-1
    endif
    n=n/2
    call imult(i2,int(p%prime,long),n)
    if (i2 .eq. 0) ovflow=.true.
  endif
end do

if (ovflow) then
  i2=0
  print "(a)", 'Integer overflow in conversion to (a/b)sqrt(c/d) form'
  call report ("Try power-of-prime representation")
end if

END SUBROUTINE rrf_to_4i

!-----------------------------------------------------------------------

SUBROUTINE char4i(k, m, m1,m2)
!  Convert K to 4-integer form (I1/I2)*sqrt(I3/I4) as characters
!  in the string M, starting at M1.
!  The pointer M1 is incremented to point at the next available
!  position in the string.

TYPE(rrf), POINTER :: k
INTEGER, INTENT(IN) ::  m2
CHARACTER(LEN=*), INTENT(INOUT) :: m
INTEGER, INTENT(INOUT) :: m1

INTEGER(long) :: i1,i2,i3,i4, modi1
LOGICAL :: brackt

call rrf_to_4i (k, i1,i2,i3,i4)
!  Overflow?
if (error) return

!  Clear buffer
m(m1:m2) = ""

!  Check for zero
if (i1 .eq. 0 .or. i3 .eq. 0) then
  call put_char(m, m1,m2, "0")
  return
endif

!  Attach sign
if (i1 .lt. 0) call put_char(m, m1,m2, minus)
modi1=abs(i1)
!  Check for anything before square root
if (modi1 .ne. 1 .or. i2 .ne. 1) then
  !  Brackets only if I2, and I3 or I4, not equal to 1
  brackt=(i2 .ne. 1 .and. (i3 .ne. 1 .or. i4 .ne. 1))
  if (brackt) call put_char(m, m1,m2, bra)
  call exp_to_char(modi1,1, m, m1,m2, -1)
  !  Put in slash if I2 not equal to 1
  if (i2.ne.1) then
    call put_char(m, m1,m2, slash)
    call exp_to_char(i2,1, m,m1,m2, -1)
    if (brackt) call put_char(m, m1,m2, ket)
  endif
  if (i3 .eq. 1 .and. i4 .eq. 1) return
  call put_char(m,m1,m2, star)
endif
if (i3 .ne. 1 .or. i4 .ne. 1) then
  m(m1:m1+4)="sqrt("
  m1=m1+5
  call exp_to_char(i3,1, m,m1,m2, -1)
  if (i4 .ne. 1) then
    call put_char(m, m1,m2, slash)
    call exp_to_char(i4,1, m, m1,m2, -1)
  endif
  call put_char(m,m1,m2, ")")
  return
endif
call put_char(m, m1,m2, digit(1))       !    Really?

END SUBROUTINE char4i

!-----------------------------------------------------------------------

SUBROUTINE write4i(k, buffer, m1,m2)

TYPE(rrf), POINTER :: k
INTEGER, INTENT(IN) ::  m2
CHARACTER(LEN=*), INTENT(INOUT) :: buffer
INTEGER, INTENT(INOUT) :: m1

INTEGER :: start

start=m1

if (k%factor == zero) then
  write (lout,"(a)") "zero"
  return
end if

call char4i(k, buffer, m1,m2) 
if (.not. error) then
  write (lout,"(a)") buffer(start:m1)
  if (k%factor > one) then
    write(lout,"(a)") "multiplied by unfactorized factor"
#ifdef ARPREC
    call mpwrite(lout,k%factor)
#else
    write (lout,"(i0)") k%factor
#endif
  end if
end if

END SUBROUTINE write4i

!-----------------------------------------------------------------------

SUBROUTINE rrf_to_char(k, buffer, m, mp)
!  Converts K to character form in array buffer, starting at position
!  M.  M is updated to point at the next available position
!  in the buffer.  The exponents only are output for the first
!  MP primes; subsequent entries are output as  . prime^exponent.

TYPE(rrf), POINTER :: k
CHARACTER(LEN=*), INTENT(INOUT) :: buffer
INTEGER, INTENT(INOUT) :: m
INTEGER, INTENT(IN) :: mp

TYPE(node), POINTER :: q=>null()
INTEGER :: i, p, buflength

buflength=len(buffer)
if (buflength < mp*3+2) call report("Buffer overflow", __LINE__, this_file)
if (iszero(k)) then
  buffer(m+1:m+1)=digit(0)
  m=m+2
  return
endif

q=>k%list

!  Sign
q%x=mod(q%x,4)
if (q%x .lt. 0) q%x=q%x+4
buffer(m:m)=plus
if (q%x .ge. 2) buffer(m:m)=minus
i=mod(q%x,2)
if (i .eq. 1) buffer(m+1:m+1)=lci
if (associated(q%next,q)) then  ! Sign factor only
  if (i .eq. 0) buffer(m+1:m+1)=digit(1)
  m=m+3
  return
end if
q=>q%next
m=m+3

!  Powers of first MP primes
p=1
do while (p .le. mp)
  if (prime(p) == q%prime) then
    call exp_to_char(int(q%x,long), 2, buffer, m, m+4, 1)
    q=>q%next
    if (q%prime < 0) return
  else
    call exp_to_char(0_long, 2, buffer, m, m+5, 1)
  endif
  p=p+1
end do

do
  if (buflength .lt. m+5) call report("Buffer overflow", __LINE__, this_file)

  call exp_to_char(q%prime, 1, buffer, m, buflength, -1)
  if (q%x .ne. 2) then
    buffer(m:m)=arrow
    if (buflength .lt. m+2) call report("Buffer overflow", __LINE__, this_file)
    m=m+1
    call exp_to_char(int(q%x,long), 2, buffer, m, buflength, -1)
  end if
  q=>q%next
  if (q%prime < 0) exit
  buffer(m+1:m+1)=dot
  m=m+3
end do

END SUBROUTINE rrf_to_char

!-----------------------------------------------------------------------

SUBROUTINE exp_to_char(n, x, m, m1, m2, jr)
!  Converts N/X (N, X integers) into characters in string M, starting at
!  M1 and finishing at M2, left or right justified according
!  as JR <0 or >0.  M1 is reset on exit to the next available
!  character of M.  M is assumed to be cleared to spaces on entry.
!  The value of X must be in the range 1 to 9 (it will normally be
!  1 or 2).

INTEGER(long), INTENT(IN) :: n
INTEGER, INTENT(IN) :: x, m2, jr
INTEGER, INTENT(INOUT) :: m1
CHARACTER(LEN=*), INTENT(INOUT) :: m

INTEGER :: i, j, ns, l
INTEGER(long) :: nn

if (m1 > len(m)) call report("Buffer overflow", __LINE__, this_file)
ns=m1
nn=n
l=m2+1
if (n < 0) then
  nn=-n
  m(m1:m1)=minus
  ns=m1+1
end if

!  Construct character string initially at end of buffer, working
!  backwards from m2.
if (x > 1) then
  if (mod(nn,int(x,long)) == 0) then
    nn=nn/int(x,long)
  else
    m(m2:m2)=digit(x)
    m(m2-1:m2-1)=slash
    l=l-2
  end if
end if

do
  j=mod(nn, 10_long)
  l=l-1
  if (l < ns) then
    do i=m1,m2
      m(i:i)=star
    end do
    m1=m2+1
    return
  end if
  m(l:l)=digit(j)
  nn=nn/10_long
  if (nn == 0) exit
end do

if (l == ns) then
  m1=m2+1
else if (jr > 0) then
  !  Right justify (just move sign)
  m(l-1:l-1)=m(m1:m1)
  m(m1:m1)=space
  m1=m2+1
else
  !  Left justify
  m(ns:ns+m2-l)=m(l:m2)
  l=max(ns,l)
  m(l:m2)=""
  m1=ns+m2-l+1
end if

END SUBROUTINE exp_to_char

!-----------------------------------------------------------------------

SUBROUTINE imult (i, n, l)
!  Multiply I by N**L, checking for overflow.  L >= 0.
!  If overflow occurs, I is set to zero.

INTEGER(KIND=long), INTENT(INOUT) :: i
INTEGER(KIND=long), INTENT(IN) :: n
INTEGER, INTENT(IN) :: l
INTEGER(KIND=long) :: j, m

if (l == 0) return
m=max_long_int/n
do j=1,l
  if (abs(i) > m) then
    i=0
    exit
  end if
  i=i*n
end do

END SUBROUTINE imult

!-----------------------------------------------------------------------

SUBROUTINE put_char(buffer, m1,m2, ch)

!  Insert character ch at position m1 in the buffer. There is an error
!  if m1 exceeds, which may but need not be the end of the buffer.

CHARACTER(LEN=*), INTENT(INOUT) :: buffer
CHARACTER(LEN=1), INTENT(IN) :: ch
INTEGER, INTENT(INOUT) :: m1
INTEGER, INTENT(IN) :: m2

if (m1 > m2) call report("Buffer overflow", __LINE__, this_file)
buffer(m1:m1)=ch
m1=m1+1

END SUBROUTINE put_char

!-----------------------------------------------------------------------

SUBROUTINE setf(n,k)
!  Check that factorial(N) exists, and calculate it if it does not,
!  where N=abs(I).  If optional argument K is present, set RRF K to
!  factorial(N).

INTEGER, INTENT(IN) :: n
TYPE(rrf), POINTER, OPTIONAL :: k
INTEGER :: m

if (n == 0) then
  if (present(k)) call unit(k)
  return
endif

if (n > maxfact) then
  call report('Reference to untabulated factorial')
  if (present(k)) call clear(k)
  return
endif

if (factorial(n)%unset) then
!  Factorial(N) is not tabulated; calculate it.
  m=n-1
  do while (factorial(m)%unset)
    m=m-1
  end do
  do while (m < n)
    m=m+1
    call copy(factorial(m-1)%rrf,factorial(m)%rrf)
    call multi(factorial(m)%rrf,int(m,long),1)
    factorial(m)%unset=.false.
  end do
!       write (lout,'(a,i3,a)') 'Factorial', N, ' tabulated'
endif
if (present(k)) call copy(factorial(n)%rrf,k)

END SUBROUTINE setf

!-----------------------------------------------------------------------

SUBROUTINE multf (k, n, mexp)
!  Multiply K by (factorial N)**mexp

TYPE(rrf), POINTER :: k
INTEGER, INTENT(IN) :: n
INTEGER, INTENT(IN), OPTIONAL ::  mexp

if (n .le. 1) return
if (n > maxfact) then
  call report('Reference to untabulated factorial')
  call clear(k)
  return
endif
!  Check that the factorial is tabulated
call setf(n)
call mult (k, factorial(n)%rrf, mexp)

END SUBROUTINE multf


END MODULE rrf_module
